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Abstract 


Recent astronomical observations show that the universe may be anisotropic on large scales. 
The Union2 Snla data hint that the universe has a preferred direction. If such a cosmological 
privileged axis indeed exists, one has to consider an anisotropic expanding Universe instead of the 
isotropic cosmological model. In this paper, we present a detailed analysis of the dark energy dipole 
in f(R , T ) = fi(R) + f 2 (T) Cosmological Model using three types of dipole fit (DF) method which 
are (I)dipole + monopole fitting for distance modulus(DMFDM), (Il)dipole + monopole fitting 
for luminosity distance(DMFLD) and (III) general dipole fitting for luminosity distance(GDFLD). 
We have found the maximum anisotropic deviation direction for (DMFDM) method as (l,b) = 
(315±||, —23±^), for (DMFLD) as {l,b) = (l,b) = (315±||, —23±i|), and for (GDFLD) method 
as (l,b) = (317±||, —23±i|) which are located very close to each other. We compare our model 
with the CPL, A.CDM and ujCDM models. Constraints on (l,b) in f(R,T ) model are not much 
different from the cases of the CPL , A CDM and uiCDM models. Moreover, the results are 
consistent with other studies. 
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1. INTRODUCTION 


Cosmological principle is one of the basic assumptions of modern cosmology. According 
to the cosmological principle, the Universe is homogenous and isotropic on scales larger 
than a few hundred Mpc, which is consistent with currently observational data sets such as 
the Cosmic Microwave Background (CMB) radiation data from the Wilkinson Microwave 
Anisotropy Probe (WMAP)( 38]- 60|). However, recent observational evidence included 
Large Scale Velocity Flows (DarkFlow) (|29j) anisotropy in the Values of the Fine Structure 
Constant a ( a Dipole) ( 61], 4l|), anisotropy in Accelerating Expansion Rate ( Dark Energy 


Dipole) ( 


Q. El, 


and other effects 


47j-j5]-46j) indicate that the LIniverse may be 


anisotropic on large scales. A number of authors have investigated the anisotropies of 


the cosmic acceleration ( 32 
several groups such as |5l] 


12]), which was motivated in several aspects. In particular, 
371 1 have applied the hemisphere comparison method to study 
the anisotropy of A CDM, ujCDM and the dark energy model with CPL parametrization. 
Some previous works payed attention to study the anisotropic expansion of the universe 
(Dark Energy Dipole) using the SNIa data and found statistically significant evidence for 
anisotropies. 


More recently, j|J have applied the hemisphere comparison method to the standard A CDM 
model and found that the hemisphere of maximum accelerating expansion is in the direction 
(/, b) = (309^23,18+n) with Union2 data. |5l] took use of the hemisphere comparison 
method to fit the A CDM model to the supernovae data on several pairs of opposite 
remispheres, and a statistically significant preferred axis was found. 


3 have investigated the anisotropic Cosmological mod 
They found the preferred direction as (/, b) = (306,-18). 


el in the Randers space-time. 
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have taken the deceleration 


oarameter q 0 as the diagnostic to quantify the anisotropy level in the uCDM model. 


111 ] constructed a direction-dependent dark energy model based on the isotropic background 


described by the A CDM, ujCDM and CPL models and employed the Union2 dataset to 
constrain the anisotropy direction and strength of modulation. They found the best-fitting 
value of the maximum deviation direction from the isotropic background is not sensitive to 
the details of isotropic dark energy models. 


591 ] have studied dipolar anisotropic expansion with cosmographic parameters. They found 


(l,b) = (309°, —8.6 C 
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62 ] chose two simple cosmological models, KCDM and ojCDM for the hemisphere com¬ 
parison approach, and KCDM for the dipole fit. In the first approach, they used the matter 
density and the equation of state of dark energy as the diagnostic qualities in the KCDM 
and ojCDM , respectively. In the second method, they employed distance modulus as the 
diagnostic quality in KCDM. They found a preferred direction of (l,b) = (307°, —14°). 

In testing for anisotropy or consistency with isotropy, we can ask which cosmological probes 
are most sensitive in what redshift ranges to such a hypothetical anisotropy, i.e. what 
constraints could be put on angular variations in the local dark energy equation of state. 
We cannot make a convincing conclusion from only one dataset, model or method about 
the origin of the anisotropy. Anisotropy may come from systematic uncertainty, as well as 
the intrinsic property of the universe. If the privileged axes derived from different datasets, 
different methods and different cosmological models are close to each other, we can safely 
conclude that anisotropy is an intrinsic property of the Universe. As we mentioned 
above several studies payed attention to find a preferred axis of the Universe in isotropic 
background described by the KCDM , ojCDM and CPL models; however, possibility 
of existence a privileged axis for the Universe may enhance if modified cosmological 
models such as f(R,T) could predict it and produce closely result. Searching for preferred 
cosmological axis of the Universe, we focus on generalized gravity model f(R,T). 


Recent 


galaxies ( 


& 


the observations of high redshift type la supernovae, the surveys of clusters of 


56]- 49]), Sloan digital sky survey (SDSS)( [l ]) and Chandra X-ray observatory 


reveal the universe accelerating expansion and that the density of matter is very much 
less than the critical density. Also, the observations of Cosmic Microwave Background 
(CMB) anisotropies (jzi) indicate that the universe is flat and the total energy density is 
very close to the critical one (53]). The observations though determines basic cosmological 
parameters with high precisions and strongly indicates that the universe presently is 
dominated by a smoothly distributed and slowly varying dark energy (DE) component, 
but at the same time they poses a serious problem about the origin of DE (j^]). The 
most The cosmological constant as the best candidate for explaining cosmic acceleration 
in literature faces serious problems such as fine-tuning and a huge discrepancy between 


theory and observations 


20]- 58]). On the other hand, modification of the geometrical 


part of the Einstein-Hilbert action by replacing an arbitrary function of the Ricci scalar R 
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45|) has constructed well-developed dark energy models. This phenomenological approach 
is called as the Modihed Gravity. Using the Modihed Gravity we can strongly explain 
the rotation curves of galaxies, the motion of galaxy clusters, the Bullet Cluster, and 
cosm olog ical observations without the use of dark matter or Einsteins cosmological constant 


44]- 30]). Cosmic inflation, mimic behavior of dark matter and current cosmic acceleration 


being compatible with the observational data are other successful predications of the f(R) 


theories 


44]-y) 


A generalization of f(R) modihed theories of gravity was proposed in jo] studies, by 
including in the theory an explicit coupling of an arbitrary function of the Ricci scalar 
R with the matter Lagrangian density L m . A specific application of the latter f(R, L m ) 


gravity was proposed in 


481 ] studies, which may be considered a relativistically covariant 


model of interacting dark energy, based on the principle of least action. The cosmological 
constant in the gravitational Lagrangian is a function of the trace of the stress-energy 
tensor, and consequently the model was denoted A(T) gravity. It was argued that recent 
cosmological data favor a variable cosmological constant, which are consistent with A(T) 
gravity, without the need to specify an exact form of the function A(T) ( 4.8]). A(T) gravity 
is more general than the Palatini f(R) gravity, and reduces to the latter when we neglect 
the pressure of the matter. 


In this paper, a class of the Modihed Gravity theories in which the gravitational action 
contains a general function f(R,T), where R denotes the Ricci scalar and T is the trace 
of the energy-moment urn tensor, has been considered. 33] introduced this type of the 
Modihed Gravity,/ (R, T), which obtained significant outcomes: the reconstruction of 
cosmological solutions, where late-time acceleration was accomplished by 34] and the 
energy conditions was analyzed by Q. (^2] studied the thermodynamics of Friedmann- 
Lematre-Robertson-Walker (FLRW) spacetimes. Moreover, the occurrence possibility of 
future singularities was studied by (j35(). Besides these achievements, a serious shortcoming 
in this kind of theory has been the non-conservation of the energy-momentum tensor. 
To circumvent this problem, in this paper, we show that f(R,T) functions can always 
be constructed in a way to be consistent with the energy-momentum tensor standard 
conservation. In this regard, we can assume separable algebraic functions of the form 
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f(R,T) = fi(R) + fiiT) in which the function / 2 (T) is obtained by imposing the conser¬ 
vation of the energy-momentum tensor. Now, in order to search for dipolar asymmetry, we 
construct an anisotropic dark energy model and aim to detect the maximum anisotropy 
direction. Furthermore, we consider the impact of redshift on the direction by using 
the redshift tomography method, with the Union2 data. Finally, we compare our results 
for the f(R,T) model with AC DM, luCDM and CPL models and also some previous studies. 


The paper is structured as follows: In Section 2, we obtain the field equations of 
f(R,T) gravity and analyze the stability of the dynamical system of the f(R,T) model. In 
section 3, we discuss free parameters of the model in some detail and constrain these free 
parameters using observational data. In section 4, we investigate the scalar perturbations 
in f(R, T) = fi(R) + / 2 (F) type theories. I 11 order to search for Dark Energy Dipole in the 
model using observational data, in section 5 we describe some important anisotropy models 
and method. Then we introduce and extend types of Dipole-Fitting method in order to 
investigate possible anisotropy from the data. We compare the results of these three types 
of DF method used for the f(R,T) model with each other in Section 6. Moreover, in order to 
explore the possible redshift dependence of the anisotropy, we have implemented a redshift 
tomography analysis in Section 7. In section 8, we have performed the anisotropy analysis 
for CPL , A CDM and luCDM models. We have applied DMFLD method to find the 
anisotropy of A CDM and luCDM and the dark energy model with CPL parametrization 
in order to make a comparison between these models and the f(R,T) model. Finally, 
in section 9, we conclude, summarize and compare the results of this work with some 
of the recent studies of 


ini, 0, 0|, searching for evidence for a preferred cosmological axis. 


2. FIELD EQUATIONS OF f(R,T) MODEL 


The action of f(R, T) gravity is of the form 


S 


\J—gd A x 


1 

167rG 


/(R,T (m) ) + L (m) + L (rad) 


( 1 ) 


where R is Ricci scalar, f(R, T is an arbitrary function of the Ricci scalar and T( m \ 
and l/ rad ) are the Lagrangian of the dust matter and radiation, g is the determinant of the 
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metric, T ^ = g^Tffl is the trace of the energy-momentum tensor and we set c = 1. By 
varying the action ([I]), with respect to the metric tensor g the field equations are obtained 
as 


f R (R,T)R lx „--f(R,T)g lx „+ 
(v D ^v^v.)/fl(-R,r) = 
8nG + f T (R, T))TW + 8vr GT^ d \ 


where V denotes the covariant derivative and 


□ = v*v„, 


d T 


(m) 


f (p T \ — df (Ri T) oylfm) 


( 2 ) 


(3) 


contracting of equation (j2j) yields 

f R {R,T)R + 3nf R (R,T)-2f(R,T) = (8 ttG + MR,T))T. (4) 

Now, in this model, we assume the perfect fluid and the spatially flat Friedmann- 
Lemaitre-Robertson-Walker (FLRW) metric 

ds 2 = — dt 2 + a 2 (t) [dx 2 + dy 2 + dz 2 ^j , (5) 

where a(t ) is the scale factor. Let us rewrite (J2]) as a standard form similar to GR, i.e. 

_ ^ 7r< ^ , T^(rad) , T^(eff)\ (fp 

fji(R, T) ' ^ ^ ' 3" ) 1 v u / 

where 

T tf ] = ~ ft(R,T)R)g rv + ( 7 ) 

(V,V«-9^ n ) /«(«■ T ) + MR. T)T ( ™'] 

Regarding the Bianchi identity, obviously in f(R,T) gravity, the above effective energy- 
momentum tensor is not conserved. Thus, by applying the conservation of the energy- 
momentum tensor of the whole matter, i.e. V /i T / i" 1 ' 1 = 0 = V M T^ ad \ the following constraint 
must be held. That is 

3 -H(t)f T (R,T) = MR.T), (8) 


7 







where H(t) is the Hubble parameter and dot denotes the derivative with respect to the time 
t. Equations (j2D and (j3J), by assuming metric ([5]), give 


3 H 2 f R (R, T) + i (f{R, T) - f R (R, T)R ) (9) 

+3 f R (R, T)H = (87rG + f T (R, + 8vr Gp^ d \ 

as the Friedmann-like equation, and 

2 f R (R, T)H + f R (R, T ) - f R (R, T)H = (10) 

- (SttG + f T (R, T)jp^ - -7T Gp (rad) , 


as the Raychaudhuri-like equation, 
gave three classes of these models 
R + 2f(T) 
fi(R) + f 2 (T) 
fi(R) + f 2 (R)f 3 (T) 


f(R,T) = { 


In this paper, we have focused on the second class f(R,T) = fi(R) + / 2 (E). Now, by 
rewriting equations ((9]) and (fTOl) for this model, one can obtain 

R f[(R) h(R) f 2 (T) 


6H 2 Hf[(R) 6 H 2 f[{R) 6 H 2 f[(R) 

8tt f' 2 (T)p 8vrGp( rad ) 


( 11 ) 


+ 


3 H 2 f[(R) 3 H 2 f[(R) 3 H 2 f[{R) 


= 1 


and 


H f[(R) f[(R) 

H 2 2 H 2 f[ 2 Hf[(R) 

4tt Gp^ /'(T)p( m ) 16 tt Gp( rad ) 

~H 2 f[(R) ~ 2 H 2 f[{R) ~ 3 H 2 f[(R) ' 

For more simplicity, we introduce a few independent new variables as 

R 


( 12 ) 


c = m) 


e = 


/ 2 CO 


_ /i(i2) 

V = 


X = 


6H 2 f[(RY 

8vrGp (rad) 


7 ? = 


OFT 2 ’ 


(13) 


^ Tf' 2 {T) 


3 H 2 f[(RY 3 H 2 f[(RY - 3 H 2 f[(R) 

where the prime denotes differentiating with respect to the argument and we have used 
i? = 6 (H + 2H 2 ) for metric ([5]). By applying the conservation equation (jHJ) for the minimal 























combination, f(R,T) = fi(R) + / 2 (T) gives 

Tf2(T) = ~lf' 2 (T) (14) 

This constraint restricts its form to a particular one, namely, 

f{R, T ) = h{R) + + c 2 (15) 


where c\ and c 2 are constants with respect to T. The conservation of the energymomentum 
tensor also leads to the case in which the variable v is a function of £, namely,!' = |. 
Therefore, these six variables will reduce to five independent variables once the constraint 
equation dSJ) is applied. 


d <_ 

dN 

dr] 

dN 

dd_ 

dN 

dN 

dx 

dN 


-1 + C(C - ti) - 3r] ~ ^ + X, 

— + r] (A + ( — 2d ), 
a 

-^- + 2d (2 -■&), 

a 

eQ + c-2tf), 

x(C-20), 


(16) 

(17) 

(18) 

(19) 

( 20 ) 


JD Cl I t D\ 

Where N represents derivatives with respect to In a and a = ; which for constant value 

of a gives 


MR) = c 1R R“ +1 + c 2R 


( 21 ) 


where Cm and C 2 r are constants. 

Note that the second order nonlinear differential equations of systems was simplified to 
first order differential equations by introducing a few new variables, ft is interesting to 
consider the behavior of systems around the equilibrium points= 0^ = 0^ = = 

= 0) us i n S jacobian stability analysis. 

The Jacobin stability of a dynamical system can be regarded as the robustness of the system 
to small perturbations of the whole trajectory. This is a very convenable way of regarding the 
resistance of limit cycles to small perturbation of trajectories, ft gives us the possibility of 
study all the evolutional paths admissible for all initial conditions ( 24]- 27]). ft is especially 
important in cosmology where there is the problem of initial conditions. Using the dynamical 
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TABLE I: The fixed points solutions of the dynamical system problem of f(R, T) = fi(R) + 


Fixed points Coordinates (£, rj, d, x) 


eigenvalue 


Pi 

(0,0,0,0 

,1) 

(5 1 

\ 2 ’ ■ L ’ 

-1,4,4) 


P‘l 

(-4,5,0, 

0,0) 

+ 

-5,-4,- 

-3, 1( “ +1 >) 

P.1 

(1,0,0,0 

,0) 

(l,2, 

c 7 4(a- 
°’ 2’ a 

*) 

Pi 

(-1,0,0, 

0,0) 

( 3 X 

,-2,-1, 

4a+l)\ 
a y 

p 5 

(-|,0,0 

, 1,0) 


-5 -3 3 
2 ’ 2 ’ 2 

8a+5\ 

’ 2a 1 

Pr 

( 4a 

2 a 2a rj 5a 2 +201+1^ 

f 5 1 

4 a—1+&4 a— 1—^4 \ 

l 6 

y l+a ’ 

(1 +a)' 2 ’ l+a’ U ’ (1+a)^ J 


’ 2a+l 

’ 2a+l J 

Pr 

(0,-1, 2, 

0,0) 

(* 

-4,-3,- 

-3a+5.i 3a+<5?, \ 

2a 5 2a 


P 8 


( 2(1—a) 1—4a (l-4a)(l+a) n ^ ( 2(5a 2 +2a-l) 8a 2 +3a-2 l-4a 2(a 2 -l) 10a 2 +3a-4 N \ 

\ l+2a ’ a(l+2a) ’ a(l+2a) ’ U ’ U J ^ a(2a+l) ’ a(2a+l) ’ a ’ a(2a+l) ’ 2a(2a+l) ) 


Pg 

/ 3a 

l+4a: 

1+4 a 

y l+a ’ 

2(l+a) 2 ’ 

2(l+a) ' 

Pi 0 

( 3a 

5+8a 

5+8a 

\ 2(l+a) 

’ 4(l+a) 2 

’ 4(l+a) 


, 0,0 


4(l+a) 2 


(o 3 _1 — 3a4jg 3a+<52 ( 

\+ 2’ 4a(a+l) ’ 4a(a+l) J 



^3 3 1 

2 ’ 2 ’ 8a(a+l) 


(—6a 2 


9a + 8i), — 8a (^ + i) (6a 2 + 9a + ■ 


systems methods one hopes to answer the question of what is the range of initial conditions 
and parameters of the system for which the subsequent evolution is compatible with current 


J = 


(2(-d -3 -C-l 

^ + V 4 + C-2d J - 2 77 

0 ^ + 4 - 4d 

a 

e o -2£ 

V X 0 -2x 


-■& 

a 


-3 

2 

0 

0 

+ C-2d 
0 


1 \ 
0 
0 
0 

C-2d/ 


( 22 ) 


<5i = V676a 4 + 700a 3 - 55a 2 - 16a 
5 2 = \/256a 4 + 160a 3 — 31a 2 — 16m 


Si = \/25a! 2 — 16ct 
S 4 = V81a 2 + 30a - 15 


(23) 
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Fig. 1: The attractor property of the dynamical system in the 3-Dimensional phase plane. 



Fig. 2: The attractor property of the dynamical system in the 3-Dimensional phase plane. 


Evaluating the Jacobian matrix at the steady state and computing the corresponding 
eigenvalues of them, we can investigate stability or instability based on the real parts of the 
eigenvalues. Table [D shows the property of critical points of dynamical system. Fig.l and 
Fig. 2 demonstrate the attractor property of the dynamical system in the 3-Dimensional 
phase plane from different perspectives. Imprecisely speaking, the trajectories of the phase 
space approach to a fixed point if all eigenvalues get negative values, and recede from a 
fixed point if all eigenvalues have positive values. However, the fixed points occurring in 
the former and the latter sets are called the stable and unstable points, respectively. The 
fixed points with both positive and negative eigenvalues are called saddle points, and those 
trajectories which approach to a saddle fixed point along some eigenvectors may recede 
from it along some other eigenvectors. 
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3. CONSTRAIN ON PARAMETERS OF THE MODEL 


In pervious section, we investigated stability of dynamical system by introducing 
the dimensionless parameters {C, ^ 7 , ^ 1 ? if is obvious that the critical points and 

eigenvalues dependent only on the free parameter a (see Table [I]). Also, we can see from 
equation [16] to [ 20 ] that parameter a is the only parameter which has been explicitly revealed 
in the set of equations and has directly affected the dynamical system; however, there 
are some parameters such as {Cm, Cm, C 2 R, C 2 T} which have not appeared in the set of 
equations [16] to [ 20 ] and have been masked by the dimensionless new variables; but they can 
affect on dynamics of the system. In fact, the new variables dependent on these parameters. 
For example, we can rewrite the variable 77 as p = ~ 6 H^f(R) = ~ (l+» + §^i?,A a+1 i j 
We can see that this variable dependents on Cm, C 2 r and a. It is important to note 
that in a dynamical system with a set of equations both free parameters and initial 
conditions determine the dynamics of the system. Free parameters affect critical points and 
initial conditions affect the trajectories of variables in phase space. Here, the parameter 
a is the only free parameter. Although the parameters {Cm, Cm, C 2 r , C 2 t} have not 
appeared in the equations explicitly, they can affect value of initial conditions. I 11 order to 
study the effect of these parameters and constrain them with observation, we reveal these 
parameters in new set of equations. I 11 this respect, we introduce some other new variables as 


X! = H : X2 = R, x 3 = p m , x 4 = P rad , *5 = MR, T),xq — MR , T) (24) 


Equation (j9j) . (fT5l) and (l 2 Tj) give the following constraints between /r(R, T) and the variables 
aqto £5, namely 



(25) 


( 26 ) 
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Fig. 3: One dimensional likelihood for parameters a,C ir, C 2 R, C\t and C 2 T 


Now,using (|25|) and ([26]) . for the autonomous equations of motions, we obtain 


dXi _ (c 1 B (a + l))° ~ 12X 1 
dN 6 x 1 

dx 3 

dN 
dx 4 

dN 


= ~3x 3 , 
= -4x 4 , 


dx 5 1 r , 

dN ~ 3x? U 


y/QCiTs Cm 

-)x 3 + x 4 - 


^5 \ Q: + 1 


- x 5 3aq - - 


12 a;i 

x 5 


2 Cmip- + 1 ) 


)5 ) - ~(C 2 R + C 2 t) - ClT ^ Xl \ 

2 K C 1 R {oi+iy ) 2 [2R 2T) 2 1 


(27) 

(28) 
(29) 


(30) 


Therefore, we have a dynamical system with fonr independent variables and five free 
parameters(a,C , i/j , C 2 r,C\t and Cit- Note that this system of equations is corresponding 
to equations (USD to flUD. We have best-fitted these parameters using SNe la data by y 2 
method. The likelihood for these parameters have been shown in Figl3] 


4. SCALAR PERTURBATIONS IN f(R, T) = /i(R) + f 2 (T) TYPE THEORIES 

Let us consider the scalar perturbations of a flat FRW metric in the longitudinal gauge 
and in conformal time: 


ds 2 = a 2 (r?) [(1 + 2<f>)d?? 2 - (1 - 2T)dx 2 ] , 


(31) 
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where $ = $(r/, x) and T = d/ ( 77 , x) are the scalar perturbations. The components of 
perturbed energy-momentum tensor in this gauge are given by 


6T° 0 = 6p = p 0 6, 5T) =-6p 6) =-c 2 sPo 6)S, ST* = -5T 0 = - (l + c 2 ) p 0 d iV , (32) 


where v denotes the potential for the velocity perturbations. The first order perturbed 
equations in a dust matter dominated universe, c 2 = 0 will be obtained as 4] 


<f> — T = , 

/ IRo 


(33) 


with 


5R = —- 


a 2 L 


3T" + 6 (V + V 2 ) <f> + 3 U ($' + 3T') - k 2 (<f> - 2 T) 


(34) 


3V ($' + T') + k 2 ($ + T) + 3^'T - (m' - 6 V 2 ) <f> f 1Ro + (9V$ - 3+ 3T') f' 1Ro 

= -a 2 5p 0 [K 2 - / 2 ToI35) 


<f>" + T" + 3V ($' + T') + 37f'<f> + (V + 2V 2 ) T / lflo + (3H4> — "HT + 3$') f[ Ro 

+ (3* - *0 /('«, = ^Po/m . (36) 


(24 - 1) + 4' + 4' + n (4 + 4) = —a 2 vp 0 (K 2 - f 2To ) , 


(37) 


- fc 2 u - 3T' = 0 (38) 

and 

foT 

$ + Vv + u' = n . ^ ° , (3Hu - 5) (39) 

2(e 2 - / 27 b) 

where r 2 = 87 tG, the prime holds for the derivative with respect to rj, V = a'/a and the 
subscript 0 holds for unperturbed background quantities: R 0 denotes the scalar curvature 
corresponding to the unperturbed metric, p 0 the unperturbed energy density, with / 10 = 
fi(Ro), fiR 0 = h/i(/?,())/d/?o, /20 = f2(T 0 ), / 2 t 0 = d/ 2 (T 0 )/dT 0 and c 2 = po/po , fiR 0 R 0 = 

d 2 /!(i? 0 )/di? 2 
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4.1. Solution of the equations using dynamical system 


The complete set of equations that describes the general linear perturbations for the 
model have been presented in pervious section. These equations are a set of nonlinear 
second order differential equations with a large number of variable for which there is no 
analytical solution except for simplest cases and only numerical analysis can be performed. 
Our purpose is to convert second order differential equation to first order by introducing 
some new variables. There are various reasons for doing this, one being that a first order 
system is much easier to solve numerically. Also, it allows us to investigate the behavior 
of the system in phase space. Phase planes are useful in visualizing the behavior of the 
system particularly in oscillatory systems where the phase paths can ’’spiral in” towards 
zero, ’’spiral out” towards infinity, or reach neutrally stable situations called centres. This 
is a useful method to determine whether dynamics of a system are stable or not. 

The structure of phase space of the field equations is simplified by defining a few variables 
and parameters. These variables are generally defined as 




Xi = 



k 

X2 = 

H 

X3 = 

flRo 



5 

X4 = 

¥ 


Poa 2 

X5 = 

flRo 


/ 2 T 0 

X6 = 

V? 


T 

X7 = 



T 

X8 = 

¥ 


(40) 

(41) 

(42) 

(43) 

(44) 

(45) 

(46) 

(47) 
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Now, for the autonomous equations of motions, we obtain 


dxi 
dN 
dx 2 
dN 
dx 3 
dN 
dx 4 
dW 

dx 5 

dN 

dN 

dxi 

dN 

dx 8 

dN 


r -* xi - xi£ 

(48) 

—ke 

(49) 

P~xl~ £Xz 

(50) 

n - xax\ 

(51) 

— X5 — X5X3 

(52) 

3 

2%e - 2exe 

(53) 

^ — X7X1 — &Xi 

(54) 

X7 - XsXi 

(55) 


Where N = Ina thus, Also, we have used the following parameters 


w 

n 2 

= £ 

(56) 

w 

= r 

(57) 

vp" 

§u 2 

_ 77 

(58) 

& 

= n 

(59) 

2a 2 /i Ro 
flRo Rq 

= n 

(60) 


After some calculation from equations ([33]) - ([39]) . we can obtain the above parameters in 
terms of the new variables as 
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TABLE II: Critical Points of the system 


Critical Points 

Ai 

A2 

A3 

A4 

A5 

A6 

A 7 

As 

PI 

8.6 

5.7 

-1.0 

-0.0 

-3.5 

0.0 

-0.1 

-0.0 

P2 

8.6 

5.7 

-1.0 

-0.1 

0.0 

0.0 

-8.8 

-1.0 

P3 

8.6 

-5.7 

-1.0 

-0.0 

-3.5 

0.0 

-0.1 

-0.0 

P4 

8.6 

-5.7 

-1.0 

-0.1 

0.0 

0.0 

-8.8 

-1.0 


r 

n 

n 


As 

2 


Ai + X7 + \xl{ 1 + As) + (3 - AS + At)A5 - -^AsA4(k 2 A2 - k Xe) 

Ai + A7 + Ja^ 1 + As) + (3 - As + Az)A5 - j^A5A4(k 2 A2 - ^Ae) 
As L 3 k 2 

Xi ~ 3X7 + ^A 2 - A7Ai + ^A 2 (! - 2 As) + ^(1 - J^xlxs) 

-S - 3e(l + ^Xs) - 3 xi - 3x7 - 2xs - (3 - Xs + 3 Xi)A5 + /?(as - 3) 


:A4A5A6 


~k 2 xl [(2 - Xs)A3 + Ai + X7 + 1 + As] 


As(« 2 a1 - k 2 A 6 ) 


+ 3 X7 


k 2 

3a(l + e) —7t 
A2 


(61) 

(62) 

(63) 

(64) 

(65) 


Where we have supposed that (3 = —^ 2 . By substituting equations (l6Tf - fl65|) into 
equations (j48l) - (l55lh the complete set of equations that describes the behavior of the system 
in terms of new variables will be provided. 

In general ,tlie critical points and eigenvalues of the system will be obtained in terms of /?, k. 
Here, we have obtained critical points of the system for /3 — l,k — 0.3(see table 071). The 
corresponding eigenvalues are as: 
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Ev 1 = 


/-11.7 + 86.6i\ 
-ll.7-86.6i 
—7.0 + 6.97 
-7.0 - 6.97 
-8.7 
2.3 

8.1 x 10“ 8 
1.5 


Ev 2 = 


-20.3 
2.8 

—9.9 + 5.97 
—9.9 — 5.97 
-2.5 
-7.3 

-3.1 x 10~ 7 
1.5 


\ 


Ev 3 = 



/-11.7 + 86.6A 
— 11.7 — 86.67 
—7.0 + 6.97 
-7.0 - 6.97 
-8.7 
2.3 

8.1 x 10~ 8 
1.5 


/ 


, Ev 4 = 


-20.3 
2.8 

—9.9 + 5.97 
-9.9 - 5.97 
-2.5 
-7.3 

-3.1 x 10" 7 
1.5 

( 66 ) 


\ 



Fig. 4: Attractor behavior of the system for /3 — 1, k — 0.3 


Due to the fact that there are complex values in some matrix elements of the eigenvalues, 
the dynamical system shows attractor behavior. The attractor behavior of the system has 
been shown in Fig. 4. Note that the attractor behavior in phase space implies that the 
system oscillates and moves toward steady state in a critical point. 
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The oscillating behavior of the system has been shown in Fig. 5. Moreover, we are 
interested in behavior of the parameters S , T, $. Therefore, we can reconstruct them from 
new variables as 


1 dd> 

1 dT 

§dN 

1 d5 

$dN 

Fig. 6. shows the oscillating behavior of 


Xi 

(67) 

X7 

(68) 


(69) 


parameters <5, <f>, T. 


4.2. Solution for^a£o ^ 0 

JIRq 


In 


this section, we solve the equations for the ^° R ° —> 0. Note that this limit in f(R) 

JIRq 


theories corresponds to the large scalaron mass which defined as 

a 2 fm 0 


23] 43] [28]; 


m 2 = 


3flR 0 R 0 


(70) 


[ref]. 

Applying this condition to equation (133|) yields T = <f>. Therefore, the equations (l33l) - (l39| 
are simplified as follows 


19 


























Fig. 6: The oscillating behavior of the parameters (<5, < h,'I / ). 


(6 VU + 2 k 2 <S> + 6^ 2 $)/ liJo + (6+ 3$')/^ 0 = -a 2 5p 0 [« 2 - f*r 0 ] (71) 

(2$" + + m'§ + 2"H 2 $) f 1Ro + {2U$ + 3$') f[ Ro + 2= ^a 2 5p 0 f 2 T o (72) 

(2$' + 2$^)/ 1Ro + = -a 2 V p 0 (« 2 - M) (73) 


5' - k 2 v - 3$' = 0 


(74) 


and 


$ + T~Lv + v' 


/2T 0 

2(k 2 - /27b) 


(37ft; 


8) 


(75) 


From equation (1711) we have 

i + 21 + 21 + TIl + J222l_ = i_«2 fzd + 

3W 2 W/ik, 2 •mnfw, 6if 1R „\H 2 H 2 ) 


(76) 


Hence, the autonomous Equation of Motion for the independent variables can be obtained 
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via 


dxi 

dN 
dx 2 
dN 

dx 3 

dN 

C?X4 

dN 

dN 

dxe 

dN 


r - Xi - exi 

(77) 

,n' 

~ k W 

(78) 

P~xl~ £ A'3 

(79) 

n - xaxi 

(80) 

X5 ~ X5Xl 

(81) 

3 

2^6 - 2exe 

(82) 


Also, from equation (ITS]) we obtain 

1 + Xi + + X3 + gXiX3 = frj^X 4Xb(k 2 X6 — k 2 x 2 ) (83) 

Applying constraint (183|) , the system reduces to a system with five independent variables 

0.01027895408205357 

0.01027895408205356 

0.01027895408205355 

0.01027895408205354 

0.01027895408205353 

0.01027895408205352 

0.01027895408205351 

13.5 14 14.5 15 15.5 16 

ln(a) 

Fig. 7: Fluctuation of the parameters <f>, 5 for f(R,T ) model when —> 0 



0.00944209178407320 


0.00944209178407319 


0.00944209178407318 



By setting xq = 0, behavior of the dynamical system in f(R) theory will be provided. 
Here, we have plotted two dimensional, three dimensional phase space and evolution of 
variables for f(R) theory in Fig. 8. 
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1.5 -1.4 -1.3 1.2 


1.1 


0.9 



Fig. 9: Attractor property and oscillating behavior of the dynamical system for f(R) theory 


when 0. 

JIRq 


Attractor property and oscillating behavior of the dynamical system in Fig. 9. shows 
that the trajectories spirals out from an unstable focus point and moves towards a steady 
state point. 

5. SEARCHING FOR DARK ENERGY DIPOLE USING OBSERVATIONAL 
DATA 

There are various ways to investigate possible anisotropy from the data. Generally 
speaking, there are three important ways: 

5.1. Modification of the Luminosity Distance Redshift Relation in a Specific 
Anisotropic Cosmological Model 

In this method, an expression is derived for the luminosity distance as a function 
of redshift in a specific anisotropic cosmological model. Many anisotropic cosmological 
models with modified luminosity distances have been proposed to match observations. 
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TABLE III: The incomplete list of modified luminosity distance for some Anisotropic Cosmological 
Models 



model 

modified luminosity distance 

Ref 

1 

scalar perturbation 

(O t —(11 C rz (1—dcos 9)dz 1 

[39] 

[59] 

K a L 1 -) Ho h /„ 0 Ad cos 0(l+x)5 

W^m0(l+ 2: ) 3 +l ^mO - ottIJI - 

V 6n 0 a L0 

2 

Anisotropic di, in the Finslerian space-time 

to _ n i _a c r z (l-dcose)- 1 ^ 

[ L [i + Z) Hoh J nmC -d : o S e y3+n J 

[15] 

3 

effect of peculiar velocities on di, 

d L ~ 'd-Wpec {Vpec V 0 bs)- H(z)d L \ 

[36b 

4 

wind scenario to the bulk flow 

(di = (1 + z) Jq (1 + dcos 9) = dj^{ 1 + dcos 9)) 

[9] 

\m 

5 

luminosity-distance 

in ellipsoidal universe 

dr(z9)- 1+2 f 1 dA (i-e 2 ) 1 ' 6 

- m J A ^ A 2 H (1 _ e 2 cos6 ,)l/2 

-i . 1 (1— e 2 sind) 1 / 2 

1 + Z = A' (l- e 2)l/3 

[12] 

6 

measured (perturbed) luminosity-distance 

Dl = (1 + 2h.v s )DQL ,u s =peculiar velocities 

[21] 

7 

Bianchi I Cosmology 

1 + z(9,4>, a, b, c) - [ ( a o) 
sin 2 9 sin 2 0 + 

iL{n) ~ to da35 

sin 2 9 cos 2 4> 

CO s 2 9 

C 

dz' 

2 , ) 3 +^de,0(l+2') 3[1+U ' (ft)1 

[55] 


Table Hill shows modified luminosity distance for some of these models as an incomplete list. 


The Bianchi I type cosmological model 


12], 50]) and the Rinders-Finsler cosmological 


model (|l4j-fl5|) are two models which are consistent with the SNe la data. A scalar 
perturbation of the KCDM model may also break the spherical symmetry of the Universe 


such that a preferred axis arises. 


131] . 1391 ]. 591 ] have presented a scalar perturbation for 


the AC DM model. Using a scalar perturbation for the FRW metric EU modification of the 
Luminosity Distance Redshift Relation in a specific anisotropic cosmological model will be 
obtain as [3jJ] 


d L {z, n) fa (xs - Xo)( 1 + *){l + v s • n 


v ■ n 


|Xs 

Xo 


(Xs - Xo)T-is 


(84) 


i r Xs 

V 2 (4> + T) 

^ J Xo 


{x-Xo)iXs-x) 


dx } • 


'Xo As Xo 

Here the luminosity distance is expressed in terms of the observed redshift z, and the direc- 
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tion to the source, where n denotes a unit spatial vector from the observer to the source. 
The notation ss denotes an approximate equality accurate up to first-order in the potentials 
<f>, T, and the peculiar velocities (in the conformal spacetime) v. The subscripts s and o 
refer to the source and the observer, respectively. The affine parameter y is given by 

Xs = / ~Jr\ dz + Xo , 

Jo H(z) 

This relation is appropriate for investigation of the peculiar velocities. It should be 
mentioned that the aim of our study is to investigate Dark Energy Dipole in the model. 


5.2. Hemisphere Comparison (HC) method 


The HC method divides the data points into two subsets according to their position 
in the sky and fit the subsets to an isotropic cosmological model (e.g.,A CDM model). 
Several groups such as Jsi|— j3Ti] have applied the hemisphere comparison method to study 


the anisotropy o 


A CDM, ujCDM and the dark energy model with CPL parametrization. 


More recently, |5j have applied the hemisphere comparison method to the standard A CDM 
model and found that the hemisphere of maximum accelerating expansion is in the direction 
of (l,b) = (309^23,18+}?). [5l| took use of the hemisphere comparison method to fit the 
A CDM model to the supernovas data on several pairs of opposite hemispheres, and a 
statistically significant preferred axis was found. Some of studies which have used this 
method are listed in Table IIVI 


5.3. Dipole-Fitting (DF) method 

Using the DF method, we can directly fit the data to a dipole (or dipole plus monopole) 
model. If the Universe is really intrinsically anisotropic and there exists a preferred 
direction, it should directly affect the expansion rate of the Universe, leading to the 
anisotropic luminosity distance and anisotropic distance modulus. In fact, this method 
corresponds to the fluctuation of the distance modulus. Anisotropic Dipole-fitting method 
has been used for searching the anisotropy of fine structure constant using quasars data on 
cosmological scale. [5) firstly applied this method to anisotropic study using SNe la dataset. 
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TABLE IV: Incomplete list of previous studies using HC and DF method 



Method 

Equation 

Ref 

1 

HemisphereComparison(HC)method 

A^Om _ 9 / ^Om,u ^0 m,d\ 

MM 

2 

HemisphereComparison(HC)method 

Ago _ r,/<?0,u— Q0,d\ 

Qo 'Qo,u+qo,d 

[10] 

3 

a(Dipole + Monopole)Fit 

47 = A cos 6 + B 

ta 

4 

Dipole + MonopoleFitting f or DistanceM odulus(D MF D M) 

(^p = dicosO + mi) 

m, 




,[59] , 




[a 

5 

GeneralizedDipoleFitting for Luminosity Distance{GDFLD) 

} g(z)(z.n ) g{z)cos6) 

m 


62] have applied this method to investigate dipolar asymmetry of the Universe. 17] have 


made a comprehensive comparison between the HC method and the DF method using the 
Union2 dataset. 

Several studies payed attention to the fluctuation of the distance modulus in order to find 
the preferred axis of the Universe using the DF method (See Table UVD . It is worth to 
mention that the anisotropic property of the Universe directly affect the luminosity distance 
and leading to the anisotropic luminosity distance. Therefore, besides the DF method 
for the distance modulus, we have used this method for Luminosity distance. We have 
explained three types of Dipole-fitting method which are based on deviation of distance 
modulus and Luminosity distance from their best values in isotropic model to investigate 
the anisotropic expansion of the LIniverse. 


A. Dipole+Monopole Fitting for Distance Modulus(DMFDM) 

The main steps of the DMFDM are shown as follows: 

(I) Convert the equatorial coordinates of SNe la to galactic coordinates 

(II) Calculate the angle of each supernova with respect to the dipole axis, which is determined 
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Fig. 10: Two dimensional likelihood for parameters (l,b) in f(R,T) model using DMFDM 
method (used y 2 analysis with 10 5 datapoints) 



= di cos 6 4- mi 


(/.*) = 


2 

\sn 

| 541 5 
I 541 

540.5 

W40 

BoS 1 ' 


Fig. 11: Union2 dataset and (1 — a) conhdence level for Dark Energy Dipole direction (l,b) in 
f(R,T) model using DMFDM method (used y 2 analysis with 10 5 datapoints) 

by 

cos di = Zi.n (85) 

where z, is the unit direction vector of the supernova, which can be expressed by using the 
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Galactic coordinate system. 


% = cos{li)sin{bi)i + sin(li)sin(bi)j + cos(bi)k (86) 

and n is the direction of dark energy dipole, which is the maximal expanding direction, 

n = cos(l)sin[b)i + sin{l)sin(b)j + cos{b)k (87) 

where (l, b ) is the Galactic coordinate direction of dipole axis 
(III) Define the angular distribution model with dipole and monopole 


= d\ cos 9i + mi 

T ) i 


( 88 ) 


where m\ and d\ denote the monopole and dipole magnitude, respectively , ft is the dis¬ 
tance modulus predicted by the isotropic f(R, T) model and p is the true luminosity distance 
of the supernova 

(IV) Fit the SNIa data by minimizing the x 2 sn value of the distance modulus. The x 2 n f or SNIa 
is obtained by comparing theoretical distance modulus with observed p obs of supernovae.we 
suppose the experiment error between each measurement is completely independent, so the 
covariance matrix can be simplified as the diagonal component, and the xl n can b e written 
as 


2 v' - SHF? 
x,n ti 

where p th {zf) is the theoretical distance modulus which it will be obtain from equation §88 j) 
as 



p th (zi) = n(zi)( 1 + m + d cos Of) 


(90) 


and p(zi) = 5log w [dL(z)] + 42.38 — 5logi 0 h also for a FRW cosmological model, one has 


d L {z) = (! + *) 


" 2 //n 


-dz!. 


(91) 


Jo H(z') 

to match the equation §UJ ]) with set of equations 1751 to \%W. we can express the equation fUT\) 
by two new differential equations as 


d(d L (z)) / - , s e 


dN 


2N 


= - d L (z ) + 


H(z) 


(92) 







Fig. 12: One dimensional likelihood for parameters(di, m\) in f(R,T) model using DMFDM 

method (used y 2 analysis with 10 5 datapoints) 



Fig. 13: Two dimensional likelihood for parameters^!,mi) in f(R,T) model using DMFDM 

method (used y 2 analysis with 10 5 datapoints) 


dN 



H{z ) (■& - 2) 


(93) 


where hq = 42.384 — 51og h , H 0 = lOO/i km.s 1 .Mpc 1 and n° bs (zi) is the measured distance 
modulus from the Union2 data. 


In this step, we employ the Union2 dataset to constrain the anisotropic dark energy 
model. The directions of the SNIa that we have used here are given in 8] work, and are 
described in the equatorial coordinates (right ascension and declination). In order to make 
comparisons with other results, we convert these coordinates to the galactic coordinates 

m (ban. 

The parameters need to be constrained are (hi, mi, l, b). Using the least y 2 n method, we can 
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find the best-fit parameters^!, mi, /, b ). The best-fit dipole direction is found to be towards 

(/, b) = (315° ± 25°, —23° ± 15°) (94) 

Fig. 10. shows the two dimensional likelihood for parameters (l, b) in f(R,T) model using 
DMFDM method. The distribution of Union2 Snla Datapoints in galactic coordinates 
along with the dark energy dipole direction (Z, b ) are shown in Fig.ll. The magnitude of the 
dipole and the monopole have obtained as 

d x = (1.4 ±0.8) x xl0" 3 ,mi = (-0.72 ±2.2) x 10" 4 (95) 


We can see that the magnitude of the monopole is one order of magnitude smaller than that 
of the dipole. This is consistent to the result of 41|, who obtained 


d x = (1.3 ±0.6) x xlO , m x = (2 ± 2.2) x 10 


->-4 


the result of 17], who fitted the data with a dipole only and obtained 


d x = (1.0 ±0.5) x xlO 


-3 


the result of ]59j for 0.015 < z < 8 with 


d x = (1.4 ± 0.6) x lO" 3 , mi = (2.7 ± 2.2) x 10 


-4 


the result of 62] with 


d x = (1.2 ± 0.5) x lO" 3 , mi = (1.9 ± 2.1) x 10 


-4 


(96) 


(97) 


(98) 


(99) 


We have obtained the likelihood function of each parameter by performing the y 2 
analysis using 10 5 data point. The results are shown in Fig.12. and Fig.13. 


B. Dipole±Monopole Fitting for Luminosity Distance (DMFLD) 

We perform a similar dipole±monopole fit using the Union2 data. Instead of () which 
corresponds to distance modulus deviations from its isotropic f(R,T) value, we use the 
luminosity distance deviation from its best fit isotropic f(R, T) value 


( A d L{z) \ 

V d L (z) ) l 


d L (z ) - d L (z) 

d L (z) 


d2cos0i ± m 2 


( 100 ) 
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Fig. 14: Two dimensional likelihood for parameters^, b) in f(R,T) model using DMFLD method 

(used x 2 analysis with 10 5 datapoints) 



Fig. 15: Union2 datapoints and (1 — a) confidence level for Dark Energy dipole direction^, b) in 
f(R,T) model using DMFLD method (used x 2 analysis with 10° datapoints) 

where, cIl(z) is the luminosity distance of the supernova in isotropic background and d^z) 
is the true luminosity distance or anisotropic luminosity distance of the supernova .therfore 
we can use the following expression 

d a r s {z) = d L (z), dt{z) = d L {z) ( 101 ) 
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Fig. 16: One dimensional likelihood for parameters (d 2 , m 2 ) in f(R,T) model using DMFLD 
method (used y 2 analysis with 10 5 datapoints) 



Fig. 17: Two dimensional likelihood for parameters (d 2 , tu 2 ) in f(R,T) model using DMFLD 
method (used y 2 analysis with 10° datapoints) 


using (1101 1) we can rewrite the equation (1 1001) as 

d a L nis (z ) = d^°{z)(d 2 cos9 l + m 2 + 1) (102) 

also n th (zi) = 5log w [d^ lls (z)\ + 42.38 — 5log w h 
we have found the best fitted dipole direction as 

(l,b) = (315° ± 37°,-23° ± 18°) (103) 

Fig.14. shows the two dimensional likelihood for parameters (l,b) in f(R,T) model using 
DMFLD method. The distribution of Union2 Snla Datapoints in galactic coordinates 
along with the dark energy dipole direction (/, b) are shown in Fig. 15. The magnitudes of 
the dipole and monopole have obtained as 

d 2 = (0.026 ± 0.014), m 2 = (-1.6 ± 5.4) x 10“ 3 (104) 
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Fig. 18: Two dimensional likelihood for parameters^, b) in f(R,T) model using GDFLD method 

(used x 2 analysis with 10 5 datapoints) 



Fig. 19: Union2 datapoints and (1 — a) confidence level for Dark Energy dipole direction^, b) in 
f(R,T) model using GDFLD method (used x 2 analysis with 10 5 datapoints) 

We have obtained the likelihood function of each parameter by performing the y 2 analysis 
using 10 5 data point. The results are shown in Fig.16. and Fig.17. 

C. Generalized Dipole Fitting for Luminosity Distance (GDFLD) 

Some studies have shown that the monopole is not significant (monopole magnitude is 
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Fig. 20: One dimensional likelihood for parameters( 50 , gi) in f(R,T) model using GDFLD 
method (used y 2 analysis with 10 5 datapoints) 



'0.003 


Fig. 21: Two dimensional likelihood for parameters (go, g\) in f(R,T) model using GDFLD 
method (used y 2 analysis with 10 5 datapoints) 

~ 10 _4 )( 59],|5]). We have also obtained (m ~ 10“ 4 ). Therefore, neglecting m and by 
considering dipole magnitude as a function of z ,the general case of Luminosity Distance 


dipole fit will be 


d L (z) - d° L (z ) 

d° L (z ) 


= g(z)cos9 = g(z)(z.n) 


(105) 


jll | hrst applied this method to AC DM model by assuming linear function 
of z as 

9{z) = 9o + 9iz: (106) 

we have found the best fitted dipole direction as 

(l,b) = (317° ±32°, -23° ±18°) (107) 


Fig. 18. shows the two dimensional likelihood for parameters (/, b) in f(R,T) model using 
GDFLD method. The distribution of Union2 Snla Datapoints in galactic coordinates along 
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Fig. 22: Comparison of (1 — a) confidence level of parameteres (l,b) for DMFDM, DMFLD 

and GDFLD models 

with the dark energy dipole direction (l,b) are shown in Fig.19. The magnitudes of the go 
and gi have obtained as 


go = (1.35 ± 1) x 1(T 3 , g± = (-0.4 ± 2) x 10" 4 (108) 

We have obtained the likelihood function of each parameter by performing the x 2 analysis 
using 10 5 data point. The results are shown in Fig.20. and Fig.21. 


6. COMPARISON OF THREE DF MODELS 

In the previous section, we have described three types of dipole-fitting (DF) method 
which has been used for statistical analysis in order to find the preferred cosmological axis 
of the Universe in f(R, T) model. These three types are as follows: 

(I) Dipole + Monopole Fitting for Distance Modulus (DMFDM), 

(II) Dipole + Monopole Fitting for Luminosity Distance (DMFLD), 

(III) General Dipole Fitting for Luminosity Distance (GDFLD). 


Several groups have applied DMFDM method to study t 
ujCDM and the dark energy model with CPL parametrization. 


le anisotropy of AC DM, 
11 ] have applied GDFLD 
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method to study the anisotropy of A CDM, ujCDM and the dark energy model with CPL 
parametrization. We have applied all of these DF methods to study privilege axis of the 
universe in f(R,T ) model. 

At first, it seems that these methods have a same origin (because of the direct relation be¬ 
tween p and cIl )• Also, the best fitted direction of preferred axis of these methods are very 
close to each other. In fact, DMFDM (( l,b ) = (315° ± 25°,—23° ± 15°)) and DMFLD 
(( l,b ) = (315° ± 37°,—23° ± 18°)) methods have resulted exactly the same value for the 
privilege axis of the universe in f(R,T) model. However, their 1 — cr confidence level are 
different. As left panel of Fig.16. shows, the (1 — a) confidence region of DMFDM is 
smaller than DMFLD (right panel of Fig.22.). Moreover, they give different values of dipole 
magnitude which is interesting to note. The dipole magnitude obtained using DMFDM 
method ( d\ = (1.4 ± 0.8) x xlO -3 ) is close to previous studies of 17], |59| |62] as it has 


been mentioned in DMFDM method section. However, the dipole magnitude obtained 
using DMFLD method (g ?2 = (0.026 ± 0.014)) is different from the value obtained using 
DMFDM method and also previous studies. Interestingly, the magnitude of anisotropy 
(^2 = (0.026 ± 0.014)) obtained using DMFLD method is approximately equal to that of 
CMB dipole. The recent released Planck data show that the dipole magnitude of CMB 
temperature fluctuations is about A = 0.070.01 (|l5|). 

There are two reasons to study the dark energy dipole of the universe using the formula 
based on deviation on Luminosity distance ( DMFLD method) instead of distance modu¬ 
lus ( DMFDM method). The first is that if dark energy has anisotropic repulsive force, it 
will directly affect the expansion rate of the Universe, leading to the anisotropic luminosity 
distance; therefore, in formulating the dipole-fitting method it is more appropriate that the 
cLl be revealed directly in the equation. The later reason is that most of the formulaes for 
modification of cIl presented in Table Hill with very small values of dipole magnitude (d -C 1) 
can be simplified as 


( 1l = d° L (1 ± dcosd), 


d. L ~ d° L 

d° L 


= ±dcosd 


(109) 


which is the same as the DF equation of DMFLD and GDFLD methods. 
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TABLE V: Constraints of the directions and amplitude of maximum anisotropy using DMFDM 
method for different redshift bins of the SNIa data. The error-bars quoted is la error. 


range 

l 

b 

mi 

di 

0 - 

0.2 

'2i;i 

10 +16 

—0 onot^r 000042 

U.UUUoD+ 0 0004 5 

0.00185;°;°°“° 

0 - 
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7. REDSHIFT TOMOGRAPHY ANALYSIS FOR THREE TYPES OF DF 
METHOD IN F(R,T) MODEL 


In order to explore the possible redshift dependence of the anisotropy, we implement 
a redshift tomography analysis, for the following redshift slices: 0-0.2, 0-0.4, 0-0.6, 0-0.8, 
0-1.0, 0-1.2, 0-1.4. the results of redshift tomography analysis for f(R,T) model using three 
types of the DF method are summarized in Table IVl IVII and IVIIl 


In order to elaborate redshift tomography figuratively, we have plotted the likelihood 
of the parameters (di,mi,li,bi), (d 2 , m 2 , l 2 , £> 2 ), {go, 9i, h, 63 ) f° r eac h redshift slice in some 
figures. The redshift tomography analyses in Fig. 23., Fig. 24. and Fig. 25. show 
that the preferred axes at different redshifts are all located in a relatively small region of 
the Galactic Hemisphere. The maximum anisotropic deviation direction is for (DMFDM) 
method as (l,b) = (135+25, 23l 34 ), for (DMFLD) as (l,b) = (135lp, 23l]J), and for 
(GDFLD) method as (l,b) = (137^2, 23l]J). Note that these directions are equivalent 
to (Z, 6 ) = (315±|i, —231^) for (DMFDM) method, (Z, b) = (315±|f, —23±i|) for (DMFLD), 
and (Z, b ) = (317; 3 2, —23l]J) f° r (GDFLD) method as the maximum axis. 
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23: The redshift tomography analysis for DMFDM method in f(R,T) model. 
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Fig. 24: The redshift tomography analysis for DMFLD method in f(R,T) model. 
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TABLE VI: Constraints of the directions and amplitude of maximum anisotropy using DMFLD 
method for different redshift bins of the SNIa data. The error-bars quoted is ler error. 
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Fig. 25: The redshift tomography analysis for GDFLD method in f(R,T) model. 


8. COMPARISON OF THE F(R,T) MODEL WITH CPL PARAMETRIZATION, 
uCDM AND kCDM MODELS 


In this section, we compare cnR a+1 + ctV—T Gravity model with CPL parametrization, 
uCDM and A CDM models. In the framework of a spatially flat Friedmann universe, the 
expansion history of the Universe is given by 

H\z) = H*[n m 0(1 + U ) 3 + (1 - n m0 )f{z)}, (110) 
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TABLE VII: Constraints of the directions and amplitude of maximum anisotropy using GDFLD 
method for different redshift bins of the SNIa data. The error-bars quoted is la error. 
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1.2 
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24.6 + ! 
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0 - 
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23 . 7 ;! 

0 00 1 20 -0 00080 
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0 OOOI -0,0015 
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Q = 


3w(z)£l x (z) + 1 


( 111 ) 


where H — ^ is the Hubble parameter, q is the deceleration parameter, f2 m o = y- is the 
current value of the normalized matter density, Sl x {z) is the normalized dark energy density 
as a function of redshift which evolves as an d 


f(z) = exp [3 f 1 ^ d A 

J n t + z 


( 112 ) 


Next, we turn to the parametrization of w(z). There are many functional forms of w(z) in 
the literature. In this work, we consider Chevallier-Polarski-Linder (CPL) parametrization 


introduced by 


18| , 40], which invokes as barotropic factor the known expression 


% 

w(z) = w 0 + w 1 —— (113) 

1 + z 

In this case, the equation of state becomes w(z = 0) = w 0 at present time and 
w(z —y oo) = Wo + W\ at earlier time. This simple parametrization is most useful if 
dark energy is important at late times and insignificant at early times. In addition to its 
simplicity, this CPL parametrization exhibits interesting properties. However, it cannot 
describe rapid variations in the equation of state. Using this functional form and Equation 
(11101) . Equation (|112p can be written analytically as 
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Fig. 26: Confidence levels for parameters (Q m o ,h) in CPL and A CDM models. 


TABLE VIII: Best fitted parameters for isotropy background 


Model 

OJ 

Mm 

Ua 

LIcpl 

U) 0 

U)l 

a 

V 2 • 

A-mm 

h 

c R R a+l + ctV^T 

- 

— 

- 

- 

- 

- 

1 + 0.01 

i - 0.01 

543.0747981 

n 7+0.0147 
u - 1 -0.0147 

ojCDM 

-1.05 

0.29 

- 

- 

- 

- 

- 

537.76250 

0.70i±8;8iH 

CPL 

- 

90+0.03 

u>zo -0.03 

- 

0.77 

-1.23 

0.18 

— 

541.0514134 

0.69641^^6 

A CDM 

- 

97 +O.O 3 

u>z/ -0.03 

0.73 

- 

- 

- 

- 

540.90726 

o.698+8:Citl 


H{z) 2 = /7o[L2 m0 (l + z) 3 + (1 — L2 m0 )(l + ~) 3(1+aj °+ wi) 
-3uj 0 z 

e M x )] 


(114) 


While in the case of utCDM model, the equation of state of dark energy is parameterized 
by a constant cu = ^; therefore, we have 


H(z ) 2 — iLg[fi m0 (l + Z Y + (1 — + z) 3 ( 1+1 Uj (US) 


Using Union2 data and by y 2 method, we have best fitted parameters Ll m 0 ,D\,h for 
A CDM and parameters L2 m0 +cpl+, <^o, cui for CPL model. Fig. 6 shows the confidence 
levels for parameters (L2 m0 , h ) in both A CDM and CPL models. For AC DM we have ob¬ 
tained L2 m o = 0.27, = 0.73 and h = 0.698. For CPL parametrization, we have obtained 
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TABLE IX: the preferred direction of the universe in A CDM , CPL and ujCDM models 


Ref 

ACDM 

ujCDM 

CPL 

this study 

307+32,-16+i| 

308;l,-i8;}6 

QnQ — 30 —17~ 16 
auy +27 , u * +17 

m 

308^22 >-16+21 

308+?®,-14“^ 

307 -21 — 1 U 32 

ou '+17> i0 +17 

[ 62 ] 

307+11, - 14+18 

307+16 > —14+10 

- 

m 

309+11, -15+11 

- 

- 


n m0 = 0.23, 12a = 0.77, h = 0.6964, u ; 0 = —1.23 and u\ = 0.14 and for ujCDM parametriza- 
tion, we have obtained f2 m o = 0.29, c j = —1.05 and h = 0.699. The results summarized in 
Tabic mm 

We have also considered the ujCDM , A CDM and the CPL parameterized dark energy mod¬ 
els as the isotropic background. We use the isotropic background dark energy parameters 
in Table IVIIII and fit our anisotropic parameters, respectively. The results are summarized 
in Table HXl Fig.27., Fig.28. and Fig.29. show the results of constraints on (l,b) which are 
not much different from the case of the f(R,T) model. This means that the best-fitting 
value of the maximum deviation direction from the isotropic background is not sensitive 
to the details of isotropic dark energy models. The best fitted trajectories of the effective 
EoS parameter in isotropic, anisotropic c^R a+1 + c+V —T gravity, CPL and ujCDM and 
A CDM models are shown in Fig. 30. Based on this, the trajectory of anisotropy is not in 
much difference from the case of the isotropy, also the best fitted trajectory of CPL and 
A CDM models are same at late time and different in the future and ujCDM , A CDM and 
CnR a+1 + CtV~T gravity have same trajectory at late time and in the future. 


9. CONCLUSIONS 

We have investigated the cosmological solutions of f(R,T ) gravity, in isotropic and 
anisotropic space-time. In both isotropic and anisotropic cases, our studies are based on 
the phasespace analysis (the dynamical system approach). In this approach, we convert 
a set of second order differential equations to a new set of first order ones by defining 
some dimensionless variables and parameters. There are various reasons for doing this: 
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Fig. 27: Union2 datapoints and (1 — a) confidence level for Dark Energy dipole direction (l,b) in 
CPL parametrization , ujCDM and A CDM models(used Monte Carlo simulation with 10 5 

datapoints) 


TABLE X: Comparison of the preferred axes of the Universe in f(R,T), CPL, uiCDM and A CDM 
models 


Model 

l 

b 

v 2 • 

A mm 

c R R a+1 + ctV^T 

317° 

1 

to 

oo 

o 

537.76250 

CPL 

309° 

-17° 

534.87158 

A CDM 

307° 

-16° 

534.88649 

uCDM 

308° 

-18° 

535.14416 


a first order system is much easier to solve numerically, and also phase planes are useful 
in visualizing the behavior of dynamical systems, especially in oscillatory systems where 
the phase paths can ’’spiral in” towards zero, and ’’spiral out” towards infinity. Moreover, 
it gives us useful information about (in)stability of the system and critical points of the 
system. 

At first, we have obtained the field equations of f(R, T) gravity in isotropy case and have 
analyzed the stability of the dynamical system for f(R , T ) = fi(R) + /^(T). Then, we have 
studied the evolution of scalar cosmological perturbations in the metric formalism. The 
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Fig. 28: Two dimensional x 2 of (l,b) for A CDM, CPL, uCDM and f(R,T) models. 

main purpose of scalar perturbations is to find explicit expressions for the parameter $, T 
and 6 in the framework of nonlinear f(R, T) model. Unfortunately, the system of equations 
for scalar perturbations is very complicated in the case of nonlinearity. It is hardly possible 
to solve it directly. Therefore, we have used phase space approach to simplify the nonlinear 
equations of the f(R, T) = fi(R) + f 2 (T) model. We have also reconstructed the parameters 
<f>, T and 8 from new variables. In the model, the evolution of matter density perturbations 
for different cases has been studied and the corresponding results have been shown in Fig. 4, 
5, 6, 7, 8 and 9. The attractor property (spiral in and out) of the system leads to an oscillat¬ 
ing behavior of the matter perturbations and other variables and parameters. This behavior 
is predictable for the critical points of the dynamical systems whose eigenvalues are complex. 

In section 5, we have investigated Dark Energy Dipole in the f(R, T) model using Dipole 
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Fig. 29: the 1 — a errors on the Dark Energy dipole direction, for A CDM, CPL, uCDM 

and f(R, T) models. 




1 + z 


Fig. 30: The best fitted trajectory of Equation of State for 


(Left) isotropic and anisotropic f(R,T) models. 


(Right) CPL, ujCDM and A CDM models. 
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Fig. 31: The direction of preferred axis in galactic coordinate. The point red 
• denotes our result, namely, (l, b) = (317.7+32,—23.7+]J). The results for 
preferred direction in other models are presented for contrast.Point x denotes 
the result of 


591 ]. point ♦ denotes the result of 


62]],point V denotes the result 
of Ht], point T denotes the result of [161 ]. 

point □ 


point ▲ denotes the result of 

n 

denotes the result of[611], 

denotes the result of 
the result of 


fl 

point o denotes the result of 


point o denotes the result of 

Q. 


Q, 


point 


411 ]. and point + denotes 


111 ]. The light green represents 
the 1 -a errors on the Dark Energy dipole direction, which includes 
the results for preferred direction in other models. 


Fitting method. There is a range of independent cosmological observations which indicate 
the existence of anisotropy axes. This appears to be one of the most likely directions which 
may lead to new fundamental physics in the coming years. These cosmological observations 
along with their preferred directions and the corresponding references are summarized in 
Table EU 

In this paper, we present a detailed analysis of the dark energy dipoles in 
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TABLE XI: Directions of Preferred axes from different cosmological observations 


Cosmological Obs. 

l 

b 

Reference 

Point 

Dark Energy dipole 

317.7° 

-23° 

This study 

red • 

Dark Energy dipole 

309.2° 

-8.6° 

[59] 

X 

Dark Energy dipole 

307.1° 

-14.3° 

[62] 

♦ 

Dark Energy dipole 

306° 

-13° 

[in 

V 

Dark Energy dipole 

314° 

-11° 

m 

T 

a dipole 

333° 

-12° 

m 

▲ 

Dark Energy dipole 

306° 

-18° 

M 


a dipole 

331° 

-14° 

[61] 

□ 

Dark Energy dipole 

304° 

-27° 

M 

o 

Dark Energy dipole 

309° 

-18° 

m 

■ 

Dark Energy dipole 

306° 

-13° 

[in 

+ 


f(R,T) = fi(R) + / 2 (F) cosmological model using three types of dipole-fitting (DF) 
method which are (I)dipole + monopole fitting for distance modulus (DMFDM), (Il)dipole 
+ monopole fitting for luminosity distance (DMFLD) and (III) general dipole fitting for 
luminosity distance (GDFLD). 

Several groups have applied DMFDM method to study the anisotropy of AC DM, 
uCDM and the dark energy model with CPL parametrization. [ll] have applied GDFLD 
method to study the anisotropy of A CDM, ujCDM and the dark energy model with CPL 
parametrization. We have applied all of these DF methods to study privilege axis of the 
universe in f(R,T) model. 

At first, it seems that these methods have a same origin (because of the direct relation 
between /i and di). Also, the best fitted direction of preferred axis of these methods are 
very close to each other. In fact, DMFDM ((/, b ) = (315° ±25°, —23° ±15°)) and DMFLD 
(( l,b ) = (315° ± 37°,—23° ± 18°)) methods have resulted exactly the same value for the 
privilege axis of the universe in f(R,T) model. However, their 1 — a confidence level are 
different (Fig.22.). The (1 — er) confidence region of DMFDM is smaller than DMFLD. 
Moreover, they give different values of dipole magnitude which is interesting to note. The 
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dipole magnitude obtained using, DMFDM method (d\ = (1.4 ± 0.8) x xlO 3 ) is close to 


previous studies of 


17], [59], 62] as it has been mentioned in DMFDM method section. 


However, the dipole magnitude obtained using DMFLD method (d 2 = (0.026 ± 0.014)) 
is different from the value obtained using DMFDM method and also previous studies. 
Interestingly, the magnitude of anisotropy (d 2 = (0.026 ± 0.014)) obtained using DMFLD 
method is approximately equal to that of CMB dipole. The recent released Planck data 
row that the dipole magnitude of CMB temperature fluctuations is about A = 0.070.01 
l^b Also, it is close to the result of ( 15]) which have obtained the magnitude of dipolar 
asymmetry as \D\ = 0.044 ± 0.018, using modified luminosity distance in anisotropic 
cosmological model in the Finsler-Randers spacetime (formula 2 of Tabic Hill) . 


Further results of this paper are as follows: 


1. We have found The maximum anisotropic deviation direction for (DMFDM) method 
as (l,b) = (315±||, —23±i|), for (DMFLD) as (l,b) = (315±||, —23±^f), and for (GDFLD) 
method as (/, b ) = (317^32, —235i§) which are located very close to each other. Also, the 
results are consistent with other studies (41], 16], 12], 19], { 7 ], [ll|). ft is interesting that 
the results of other studies are in (1 — a) confidence level of our study (see Fig. 31). 


2. The dipole directions at high and low redshifts are in agreement. (This is confirmed 
in the redshift tomography analyses, shown in Fig. 23, 24 and 25.) 


3.We have also applied DMFLD method to find the preferred direction of the universe 
in A CDM , CPL and loCDM models using y 2 method. Our results are very close to 
pervious works which studied anisotropy in these models (see Tabl dIXD . It is interesting 
that the results of constraints on (l, b ) in f(R, T) model are not much different from the 
cases of the AC DM, loCDM and CPL models (see Table IXl and Fig. 29). This means that 
the best-fitting value of the maximum deviation direction from the isotropic background is 
not sensitive to the details of isotropic dark energy models. 
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